##########################
###############Summary by subwatershed and HUC 
##############################################


####Set working directory (example:  setwd("K:/My Drive/EPA datasets"))
setwd("K:/")
setwd("Q:/")
setwd("C:/Users/jwmille7/Dropbox/EPA_WQ_Grant_Group_Files/Water Quality Data and Metadata/Final deliverables_ 07_2021/EPA datasets")

library(reshape2) 
library(dplyr)
library(lme4)
library(ggplot2)
library(ggthemes)

#Clear all variables
rm(list=ls())

#######Get the predictor data for each subwatershed
#######################################################

land = read.csv("Prediction models/Landcover2019-12-09.csv",header=TRUE,na.strings = "-9999")
land$Chickens=land$Chickens_Broilers+land$Chickens_Layers ##Join chickens
names(land)
land$prec_mm=as.numeric(as.character(land$prec_mm))
land=land[,c(2:6,10:12,16:18,20:21,24)]

#####ADjust HL for Subwatershed #7 from 2.3 to 8.8 (NEED TO DO IN MAIN CODE!!)
lc=ifelse(land$Siteno=="A7",8.8,land$res_hydload)
land$res_hydload=lc

###GEt land covers and livestock for ungauged areas
ungauge=read.csv("Prediction models/NWALT_allyears_formodel_ungauged_withHighDens_0925_2019-10-07.csv",header=TRUE,na.strings = "-9999")
ungauge_match=read.csv("Prediction models/ungauged_sitematch.csv",header=TRUE,na.strings = "-9999")
ungauge$basin=ungauge_match$Main_Watershed[match(ungauge$SiteID,ungauge_match$SiteID_NWALT2)]
land$basin=ungauge_match$Main_Watershed[match(land$Siteno,ungauge_match$SiteID_NWALT2)]

###GEt livestock for ungauged areas
ungauge.lv=read.csv("Prediction models/Ungauged_livestock2019-10-15.csv",header=TRUE,na.strings = "-9999")
ungauge.lv=ungauge.lv[,c(2:7)]
############Broiler numbers are not right for ungauged areas as no registered chicken farms from Chatham are there.

ungauge=merge(ungauge,ungauge.lv,by=c("SiteID","Year"),all.x=T,all.y=F)

####Call Layers, all chickens..... get rid of broilers... wrong for ungauged areas.
ungauge=ungauge[,c(1:2,4:18)];     ###Get rid of "X" and Broilers
colnames(ungauge)[17]="Chickens"  ##Change name

names(land)

##########Convert ungauged land covers to km2 and join as Ag, DEv, Fst 

####Match site code and drainage areas of each HUC
ungauge$da=ungauge_match$Area[match(ungauge$SiteID,ungauge_match$SiteID)]

ungauge$ag=ungauge$Pasture+ungauge$Crop
ungauge$dev=ungauge$Residential.Development+ungauge$Transportation+ungauge$Commercial.Dev
ungauge$fst=ungauge$Low.Use+ungauge$Semi.Developed+ungauge$Wetlands

####Get land covers as km2 instead of as a percent (upper case)
ungauge$Dev=(ungauge$dev/100)*ungauge$da
ungauge$Ag=(ungauge$ag/100)*ungauge$da
ungauge$Fst=(ungauge$fst/100)*ungauge$da

###Get pre and post 1980 land covers
##################################################
######GEt older nwalt values (Dev before 1980, before 2000)
ungauge.80=ungauge[ungauge$Year==1980,] ##Get 1980 values for all watersheds
ungauge.00=ungauge[ungauge$Year==2000,]  ##Get 2000 values for all watersheds

##Attach pre-1980 values, and pre-2000 values to ungauge.yr
ungauge$pre80Dev=ungauge.80$Dev[match(ungauge$SiteID,ungauge.80$SiteID)]
ungauge$pre00Dev=ungauge.00$Dev[match(ungauge$SiteID,ungauge.00$SiteID)]

####Manipulate to get pre-1980, 1980-2000, post- 2000 land covers
ungauge$Dev_pre1980=ifelse(ungauge$pre80Dev>ungauge$Dev,ungauge$Dev,ungauge$pre80Dev)
ungauge$Dev_post2000=ifelse(ungauge$pre00Dev<ungauge$Dev,ungauge$Dev-ungauge$pre00Dev,0)
ungauge$Dev_80_2000=ifelse(ungauge$pre80Dev<ungauge$Dev,ungauge$Dev-ungauge$pre80Dev-ungauge$Dev_post2000,0)


####Add precipitation data
ungauge.prec=read.csv("Prediction models/Precipitation_modelinput_ungauged_yearly_2019-10-14.csv",header=TRUE,na.strings = "-9999")
ungauge.prec=ungauge.prec[,c(2:4)];colnames(ungauge.prec)[2]="Year"

ungauge=merge(ungauge,ungauge.prec,by=c("SiteID","Year"),all.x=T,all.y=F);
colnames(ungauge)[30]="prec_mm"

###Get rid of years without precipitation
ungauge=ungauge[complete.cases(ungauge), ]

###Constant rain for all ungauged portions
#ungauge$prec_mm=1000         #####Need to get the correct prec. for ungauged areas

####Add stream retention and hydraulic loading for each ungagued watershed
ungauge$stream_rt=ungauge_match$RT[match(ungauge$SiteID,ungauge_match$SiteID_NWALT2)]
ungauge$res_hydload=ungauge_match$Hyd_load[match(ungauge$SiteID,ungauge_match$SiteID_NWALT2)]


names(ungauge)
ungauge2=ungauge[,c(1,2,22:24,27:29,30,31:32,15:17,14)]
colnames(ungauge2)[1]="Siteno"
names(ungauge2)
names(land)

total=rbind(land,ungauge2)
colnames(total)[c(3:5)]=c("Urban","Agriculture","Undeveloped")
colnames(total)[c(6:8)]=c("Urban pre-1980","Urban post-2000","Urban_80_2000")

total$prec_mm=scale(total$prec_mm) ####Normalize precipitation
summary(total$prec_mm)

######Rename basins in total dataset
#levels(total$basin)[levels(total$basin)=="HR"] <- "Haw River, Jordan Lake"
#levels(total$basin)[levels(total$basin)=="NHC"] <- "New Hope Creek, Jordan Lake"
#levels(total$basin)[levels(total$basin)=="FL"] <- "Falls Lake"

#total$basin <- factor(total$basin, levels = c("Haw River, Jordan Lake", "New Hope Creek, Jordan Lake", "Falls Lake"))###Set the order of factors

###################
###########Land cover and livestock graphs
###################
land=total[,c(15,2,3:5,6:8,12:14)]

#########Divide Chickens by 10^4 adn Hogs and Cows by 10^3
land$Chickens=land$Chickens/100000
land$Hogs=land$Hogs/1000
land$Cows=land$Cows/1000
summary(land$Chickens)
summary(land$Hogs)
summary(land$Cows)

##Aggregate by basin
####Join all HUCS in the same basin 
land.ag=land %>%
  group_by(basin,Year) %>%
  summarise_each(funs(sum(., na.rm = TRUE)))
land.ag=as.data.frame(land.ag)

land.ag$Urban_post1980=land.ag$`Urban post-2000` +land.ag$Urban_80_2000
colnames(land.ag)[12]="Urban post-1980"

melt.land <- melt(land.ag, id.vars = c("basin","Year"))
melt.land <- melt.land[melt.land[,2] %in% 1982:2017,]  ####This limits land to only 1994-2012 values
str(melt.land)

###Land covers only
###########################
#lc=melt.land[melt.land$variable =="Urban"|melt.land$variable =="Undeveloped"|melt.land$variable =="Agriculture",]


t=Sys.Date()
#png(filename=paste("Q:/Shared drives/Jordan Lake Modeling/WatershedModeling/Manuscript/Figures/Bar charts/Land covers_3basins_no retention",t,".png",sep=""), 
   # units="cm",width=18,height=16,res=300)

melt.land %>% filter(variable %in% c("Urban", "Agriculture", "Undeveloped")) %>% 
  #ggplot(., aes(x = Year, y = value,linetype=variable)) +  ##Line types
  ggplot(., aes(x = Year, y = value,colour=variable)) +   ###Color
  geom_line()+
  #geom_bar(colour="black",stat = "identity")+
  #theme_minimal()+
  theme_few()+
  theme(panel.grid.major.x = element_blank(),panel.grid.minor.x = element_blank())+
  #theme_classic()+
  facet_wrap(~basin,scales="free",ncol=1)+
  scale_linetype_manual(values=c("solid","twodash", "dotted","threedash"))+
  #scale_colour_manual(values = c("khaki2", "#e66101","#fdb863","#b2abd2","#5e3c99"))+
  theme(legend.title = element_blank())+    #scale_x_continuous(breaks=seq(0,82,1))+
  ggtitle("Land covers by basin")+
  #ylab("TN x 10 (kg/yr)")+xlab("")+
  ylab(bquote(''* ~ km^2))+xlab("")+
  theme(legend.title = element_blank(),axis.text.x = element_text(size=10),)+
  scale_x_continuous(breaks=c(1980,1985,1990,1995,2000,2005,2010,2015))

dev.off()

###Land covers with pre and post Development
###########################

t=Sys.Date()
#png(filename=paste("Q:/Shared drives/Jordan Lake Modeling/WatershedModeling/Manuscript/Figures/Bar charts/Dev_3basins_prepost",t,".png",sep=""), 
   # units="cm",width=18,height=16,res=300)

melt.land %>% filter(variable %in% c("Urban post-1980","Urban pre-1980","Agriculture","Undeveloped")) %>% 
  ggplot(., aes(x = Year, y = value,linetype=variable)) +  ##Line types
  #ggplot(., aes(x = Year, y = value,colour=variable)) +   ###Color
  geom_line()+
  #geom_bar(colour="black",stat = "identity")+
  #theme_minimal()+
  theme_few()+
  theme(panel.grid.major.x = element_blank(),panel.grid.minor.x = element_blank())+
  #theme_classic()+
  facet_wrap(~basin,scales="free",ncol=1)+
  scale_linetype_manual(values=c("longdash","dotdash","dotted","solid"))+
  #scale_fill_manual(values = c("khaki2", "#e66101","#fdb863","#b2abd2","#5e3c99"))+
  theme(legend.title = element_blank())+    #scale_x_continuous(breaks=seq(0,82,1))+
  ggtitle("Land covers by basin")+
  #ylab("TN x 10 (kg/yr)")+xlab("")+
  ylab(bquote(''* ~ km^2))+xlab("")+
  theme(legend.title = element_blank(),axis.text.x = element_text(size=10),)+
  scale_x_continuous(breaks=c(1980,1985,1990,1995,2000,2005,2010,2015))

dev.off()

######To check total values 
lc=melt.land[melt.land$variable =="Urban"|melt.land$variable =="Undeveloped"|melt.land$variable =="Agriculture",]

land.basin=dcast(lc,basin ~ Year,fun.aggregate=sum)


lc.prepost=melt.land[melt.land$variable =="Urban pre-1980"|melt.land$variable =="Urban post-1980"|melt.land$variable =="Undeveloped"|melt.land$variable =="Agriculture",]
land.basin.prepost=dcast(lc.prepost,basin~ Year,fun.aggregate=sum)


###Livestock only
###########################


#png(filename=paste("Q:/Shared drives/Jordan Lake Modeling/WatershedModeling/Manuscript/Figures/Bar charts/Livestock_basins_2",t,".png",sep=""), 
   # units="cm",width=18,height=16,res=300)
melt.land.live=melt.land %>% filter(variable %in% c("Chickens","Hogs","Cows"))

melt.land.live %>% filter(Year >1993) %>% 
  #ggplot(., aes(x = Year, y = value,linetype=variable)) +  ##Line types
  #ggplot(., aes(x = Year, y = value,linetype=basin)) +  ##Line types
  ggplot(., aes(x = Year, y = value,colour=variable)) +   ###Color
  geom_line()+
  #geom_bar(colour="black",stat = "identity")+
  #theme_minimal()+
  theme_few()+
  theme(panel.grid.major.x = element_blank(),panel.grid.minor.x = element_blank())+
  #theme_classic()+
  facet_wrap(~basin,scales="free",ncol=1)+
  #facet_wrap(~variable,scales="free",ncol=1)+
  #scale_linetype_manual(values=c("longdash","dotted","solid"))+
  #scale_linetype_manual(values=c("longdash","solid","dotted"))+
  #scale_fill_manual(values = c("khaki2", "#e66101","#fdb863","#b2abd2","#5e3c99"))+
  theme(legend.title = element_blank())+    #scale_x_continuous(breaks=seq(0,82,1))+
  #ggtitle("Livestock counts by basin")+
  #labs(x="",y=expression(paste("Chickens x ", 10^{4},"\n","Hogs and Cows x ",10^{3})))+  ##labs("Developed; pre-1980 (kg/(km2*yr)") +
  
  ylab(bquote(atop(paste("Chickens x ", 10^{4}),
                   paste("Hogs and Cows x ",10^{3}))))+
  #x = expression(paste("Volume ", m^{3},)), y = expression(paste("Volume ", m^{3}))
  #xlab(expression(paste("line1 \n line2 a" [b])))
  #ylab("TN x 10 (kg/yr)")+xlab("")+
  #ylab("count")+xlab("")+
  theme(legend.title = element_blank(),text = element_text(size=11))+
  scale_x_continuous(breaks=c(1980,1985,1990,1995,2000,2005,2010,2015))

dev.off()

#############################################
######################Add predictor data to the total values
#############################################
disch = read.csv("Prediction models/WWTP_TNloads_kg_year.csv",header=TRUE,na.strings = "-9999")
disch=disch[,c(2,4:(ncol(disch)-3))]
colnames(disch)<- c("PERMIT",seq(1994,2018,by=1))

melt.disch <- melt(disch, id.vars = c("PERMIT"))
melt.disch <- melt.disch[melt.disch[,2] %in% 1981:2017,]  ####This limits land to only 1994-2012 values
colnames(melt.disch)[2:3]<-c("Year","WWTP_TNdisch")
str(melt.disch)

disch_match = read.csv("Prediction models/Discharger_subwsd_match.csv",header=TRUE,na.strings = "-9999")

melt.disch$subwsd=disch_match$subwsd[match(melt.disch$PERMIT,disch_match$PERMIT)]
melt.disch=melt.disch[complete.cases(melt.disch), ]

melt.disch$basin=disch_match$basin[match(melt.disch$PERMIT,disch_match$PERMIT)]

#####Get graph of discharger trends through time by basin
disch.summ=melt.disch[,c(2,5,3)]

####Join all HUCS in the same basin 
disch2=disch.summ %>%
  group_by(basin,Year) %>%
  summarise_each(funs(sum(., na.rm = TRUE)))
disch2=as.data.frame(disch2)

disch2$Year=as.numeric(as.character(disch2$Year))

disch2$basin=as.factor(as.character(disch2$basin))

######Rename basins in discharger dataset
#levels(disch2$basin)[levels(disch2$basin)=="HR"] <- "Haw River, Jordan Lake"
#levels(disch2$basin)[levels(disch2$basin)=="NHC"] <- "New Hope Creek, Jordan Lake"
#levels(disch2$basin)[levels(disch2$basin)=="FL"] <- "Falls Lake"

#disch2$basin <- factor(disch2$basin, levels = c("Haw River, Jordan Lake", "New Hope Creek, Jordan Lake", "Falls Lake"))###Set the order of factors

###################
###########Discharger graph
###################

t=Sys.Date()
#png(filename=paste("Q:/Shared drives/Jordan Lake Modeling/WatershedModeling/Manuscript/Figures/Bar charts/Dischargers_3basins_no retention",t,".png",sep=""), 
    units="cm",width=18,height=16,res=300)

disch2 %>% filter(Year > 1994 & Year < 2018) %>% 
  ggplot(., aes(x = Year, y = WWTP_TNdisch / 100000)) + 
  geom_bar(colour="black",stat = "identity")+
  #theme_minimal()+
  theme_bw()+
  theme(panel.grid.major.x = element_blank(),panel.grid.minor.x = element_blank())+
  #theme_classic()+
  facet_wrap(~basin,scales="free",ncol=1)+
  #scale_fill_manual(values = c("khaki2", "#e66101","#fdb863","#b2abd2","#5e3c99"))+
  theme(legend.title = element_blank())+    #scale_x_continuous(breaks=seq(0,82,1))+
  ggtitle("TN point source discharges by basin")+
  #ylab("TN x 10 (kg/yr)")+xlab("")+
  ylab(bquote('TN x '* ~ 10^5*' (kg/yr)'))+xlab("")+
  theme(legend.title = element_blank(),axis.text.x = element_text(size=10),)+
  scale_x_continuous(breaks=c(1995,2000,2005,2010,2015))

dev.off()

###################################
#####Aggregate Dischargers by Siteno with/without retention to join to "total" 
###################################
disch.summ.site=melt.disch[,c(4,2,3)]

####Discharges without stream retention
disch.ag.site=disch.summ.site %>%
  group_by(subwsd,Year) %>%
  summarise_each(funs(sum(., na.rm = TRUE)))
disch.ag.site=as.data.frame(disch.ag.site)


##############################
############Aggregate land covers, livestock and dischargers for one figure
##############################

####Get all land covers
melt.land.lc=melt.land %>% filter(variable %in% c("Urban post-1980","Urban pre-1980","Agriculture","Undeveloped"))
melt.land.lc$source="Land use"
names(melt.land.lc)

####Get all livestock
melt.land.live=melt.land %>% filter(variable %in% c("Chickens","Hogs","Cows"))
melt.land.live$source="Livestock"

###Get all dischargers

disch2$variable="Discharger"
disch2$source="Discharger"
melt.disch=disch2[,c(1:2,4,3,5)];colnames(melt.disch)[4]="value"
melt.disch$value=melt.disch$value/100000

#####Merge all together and graph
#####################

graph3=rbind(melt.land.lc,melt.land.live,melt.disch)

graph3$source <- factor(graph3$source, levels = c("Land use", "Discharger","Livestock"))###Set the order of factors
graph3$variable <- factor(graph3$variable, levels = c("Agriculture", "Urban pre-1980","Urban post-1980","Undeveloped","Discharger","Chickens","Hogs","Cows"))###Set the order of factors


#####Add column specific to each line to control color and line type

#graph3$linetype=interaction(graph3$variable,graph3$source)
#graph3$linetype <- factor(graph3$variable, levels = c("Agriculture.Land use", "Urban pre-1980.Land use","Urban post-1980.Land use","Undeveloped.Land use","Discharger.Discharger","Chickens.Livestock","Hogs.Livestock","Cows.Livestock"))###Set the order of factors
#graph3$Year=graph3$Year-1900

t=Sys.Date()
png(filename=paste("Q:/Shared drives/Jordan Lake Modeling/WatershedModeling/Manuscript/Figures/Bar charts/All_summary_data_3basins_3sources",t,".png",sep=""), 
    units="cm",width=17,height=14,res=300)


graph3 %>% filter(Year > 1994 & Year < 2018) %>% 
  ggplot(., aes(x = Year, y = value,colour=variable, linetype=variable)) +   ###Color
  scale_linetype_manual(values=c("solid","dashed","dotdash","dotted","solid","solid","dashed","dotted"))+
  geom_line(size=1)+
  #theme_minimal()+
  theme_bw()+
  theme(panel.grid.major.x = element_blank(),panel.grid.minor.x = element_blank())+
  #theme_classic()+
  facet_grid(source ~ basin,scales="free")+
  scale_colour_manual(values = c("#fdb863","#b2abd2","#e66101","#5e3c99","black","#3182bd","#9ecae1","#d7191c"))+
  #theme(legend.title = element_blank())+    #scale_x_continuous(breaks=seq(0,82,1))+
  #ggtitle("TN point source discharges by basin")+
  #ylab(bquote(''* ~ km^2)) +
  #ylab("TN x 10 (kg/yr)")+xlab("")+
  #ylab(bquote('TN x '* ~ 10^5*' (kg/yr)'))+xlab("")+ 
  ylab(bquote(atop(paste("Chickens x ", 10^{5},"                                                                                   "),paste("Hogs and Cows x ",10^{3},"             TN x ", 10^{5}," (kg/yr)","                            ",km^{2},"              ")))) +
  xlab("")+
  #ylab(bquote('TN x '* ~ 10^5*' (kg/yr)'))+xlab("")+ 
  theme(legend.title = element_blank(),text = element_text(size=10),
         panel.spacing = unit(1.1, "lines"))+
  scale_x_continuous(breaks=c(1995,2005,2015))

dev.off()


###################################
#####Aggregate Dischargers by basin with/without retention to join for bar charts
###################################

####Discharges without stream retention
disch.ag=disch.summ %>%
  group_by(basin,Year) %>%
  summarise_each(funs(sum(., na.rm = TRUE)))
disch.ag=as.data.frame(disch.ag)

####Discharges with stream retention (Add RT and HL first)


####################################################
#################Get the model and mean posterior values to make predictions
####################################################
library("rstan") # observe startup messages

#model=readRDS("/My Drive/Post-Doc 2019/Jordan Lake Watershed Model/STAN models/Posteriors/Plots/2019-10-21_TN_prepost_1980_20000/R_data2019-10-21.rds")
####PIC_Q_ model
#model=readRDS("/My Drive/Post-Doc 2019/Jordan Lake Watershed Model/STAN models/Posteriors/Plots/2019-11-19_TN_PIC_Q_prepost_10000_1chains_GOODPIC/R_data2019-11-19.rds")
#####PIC_Q and power PIC land cover model
setwd("Q:/")
###Dec 14th, 20000
model=readRDS("/My Drive/Post-Doc 2019/Jordan Lake Watershed Model/STAN models/Posteriors/Plots/2019-12-05_TN_PIC_power_20000_3chains_Dec_14thpaper/R_data2019-12-05.rds")
##NOv. 25th, 20000
model2=readRDS("/My Drive/Post-Doc 2019/Jordan Lake Watershed Model/STAN models/Posteriors/Plots/2019-11-25_TN_PIC_power_20000_3chains_Nov29thpaper/R_data2019-11-25.rds")
###Dec 5th, 5000
model2=readRDS("/My Drive/Post-Doc 2019/Jordan Lake Watershed Model/STAN models/Posteriors/Plots/2019-12-04_TN_PIC_power_5000_3chains_Dec_14thpaper/R_data2019-12-04.rds")

#####Extract mean parameters from the model
theta = extract(model)
theta=as.data.frame(theta)
names(theta)
par.means=as.data.frame(apply(theta,2,mean))

theta2 = extract(model2)
theta2=as.data.frame(theta2)
names(theta2)
par.means2=as.data.frame(apply(theta2,2,mean))

########Get Export Coefficients
#Be_a=par.means["Be_a",];Be_d=par.means["Be_d",];Be_w=par.means["Be_w",]
Be_a=par.means["Be_a",];Be_devpre=par.means["Be_d_pre",];Be_devpost=par.means["Be_d_post",];Be_w=par.means["Be_w",]
Be_ch=par.means["Be_ch",];Be_h=par.means["Be_h",];Be_cw=par.means["Be_cw",]

########Get Precipitation Impact Coefficients
#pic_a=par.means["pic.1",];pic_d=par.means["pic.2",];pic_w=par.means["pic.3",]
pic_a=par.means["pic.1",];pic_d_pre=par.means["pic.2",];pic_d_post=par.means["pic.3",];pic_w=par.means["pic.4",]
pic_ch=par.means["pic.5",];pic_h=par.means["pic.6",];pic_cw=par.means["pic.7",]
PIC_q=par.means["PIC_q",]; #PIC_q=.07
pic_a2=par.means["pic_p.1",];pic_d_pre2=par.means["pic_p.2",];pic_d_post2=par.means["pic_p.3",];pic_w2=par.means["pic_p.4",]
#pic_ch2=par.means["pic_p.5",];pic_h2=par.means["pic_p.6",];pic_cw2=par.means["pic_p.7",]
pic_ch2=par.means["pic_p.5",];pic_h2=par.means["pic_p.6",];pic_cw2=par.means["pic_p.7",]

############Get stream retention values and Discharger Export coefficient
Sn=par.means["Sn",];Sn2=par.means["Sn2",]; Disch.dc=par.means["Be_dch",]
#Sn=.5 To check a given stream retention


#########################
#######Get expected Source TN for each subwatershed  at will reach the reservoirs
############################

######Get different values based on PIC_Q, PIC_Q2, and PIC_Q_power
#mdl=1; ###PIC_Q and PIC_Q2
mdl=2; ###PIC_Q_power

###Get precipitation for power relationship
#####Convert prec_mm to prec/mean(prec)
mean=1179.6;stdev=171.8

prec_mm2=total$prec_mm
prec_mm2=prec_mm2*stdev+mean
prec_mm_in=prec_mm2/mean(prec_mm2,na.rm=T)
summary(prec_mm_in)
quantile(prec_mm_in,c(.1,.9),na.rm=T)

############################################
#####Loop to get total contribution using either linear or power PICs for land covers
############################################

if (mdl== 1){
total$A=  (Be_a*(1+pic_a*total$prec_mm) * total$Agriculture);
#total$D = Be_d*(1+pic_d*total$prec_mm) * total$Urban;  ## Normal
total$D_pre = Be_devpre*(1+pic_d_pre*total$prec_mm) * total$`Urban pre-1980` 
total$D_post = Be_devpost*(1+pic_d_post*total$prec_mm) * (total$Urban_80_2000+total$`Urban post-2000`);  ##Pre and post
total$W = Be_w*(1+pic_w*total$prec_mm) * total$Undeveloped;

total$Ch=Be_ch*(1+pic_ch*total$prec_mm) * total$Chickens;
total$H=Be_ch*(1+pic_h*total$prec_mm) * total$Hogs;
total$Cw=Be_cw*(1+pic_cw*total$prec_mm) * total$Cows;
total$L=total$Ch+total$H+total$Cw
  }  else  {
  total$A=  (Be_a*prec_mm_in^pic_a2 * total$Agriculture);
  #total$D = Be_d*(1+pic_d*total$prec_mm) * total$Urban;  ## Normal
  total$D_pre = Be_devpre*prec_mm_in^pic_d_pre2 * total$`Urban pre-1980` 
  total$D_post = (Be_devpost*prec_mm_in^pic_d_post2) * (total$Urban_80_2000+total$`Urban post-2000`);  ##Pre and post
  total$W = Be_w*prec_mm_in^pic_w2 * total$Undeveloped;
  
  total$Ch=Be_ch*prec_mm_in^pic_ch2 * total$Chickens;
  total$H=Be_ch*prec_mm_in^pic_h2 * total$Hogs;
  total$Cw=Be_cw*prec_mm_in^pic_cw2 * total$Cows;
  total$L=total$Ch+total$H+total$Cw
    }

######ADD dishcargers aggregated by Year and Subwatershed
colnames(disch.ag.site)[1]="Siteno"
total=merge(total,disch.ag.site,by=c("Siteno","Year"),all.x=T,all.y=F);
colnames(total)[length(total)]="Dischargers"
total[,"Dischargers"][is.na(total[,"Dischargers"])] <- 0 ##Gets rid of NAs in this column for Dischargers


####################################
############Get predictions for normal, high, low precipitation values
######################################
###############This is done for all years. SHould be reported for just one year

precipitation_TN=function(input,prec){
input=total
  input$A.=  (Be_a*(1+pic_a*total$prec_mm) * total$Agriculture);
input$D = Be_d*(1+pic_d*total$prec_mm) * total$Urban;  ## Normal
#D = Be_d*(1+Bp[2]*av_prec) * D.lc.pre + Be_d*(1+Bp[4]*av_prec) * D.lc.post;  ##Pre and post
input$W = Be_w*(1+pic_w*total$prec_mm) * total$Undeveloped;

input$Ch=Be_ch*(1+pic_ch*total$prec_mm) * total$Chickens;
input$H=Be_ch*(1+pic_h*total$prec_mm) * total$Hogs;
input$Cw=Be_cw*(1+pic_cw*total$prec_mm) * total$Cows;
input$L=total$Ch+total$H+total$Cw

}

precipitation_TN_prepost=function(input,prec){
  input=total
  input$A.=  (Be_a*(1+pic_a*total$prec_mm) * total$Agriculture);
  #input$D = Be_d*(1+pic_d*total$prec_mm) * total$Urban;  ## Normal
  input$D_pre = Be_devpre*(1+pic_d_pre*total$prec_mm) * total$`Urban pre-1980`;  ## Normal
  input$D_post = Be_devpost*(1+pic_d_post*total$prec_mm) * (total$Urban_80_2000+total$`Urban post-2000`);  ##Pre and post
  #D = Be_d*(1+Bp[2]*av_prec) * D.lc.pre + Be_d*(1+Bp[4]*av_prec) * D.lc.post;  ##Pre and post
  input$W = Be_w*(1+pic_w*total$prec_mm) * total$Undeveloped;
  
  input$Ch=Be_ch*(1+pic_ch*total$prec_mm) * total$Chickens;
  input$H=Be_ch*(1+pic_h*total$prec_mm) * total$Hogs;
  input$Cw=Be_cw*(1+pic_cw*total$prec_mm) * total$Cows;
  input$L=total$Ch+total$H+total$Cw
  
}

precipitation_TN_prepost_power=function(input,prec){
  input=total
  input$A.=  (Be_a*prec_mm_in^pic_a2 * total$Agriculture);
  #input$D = Be_d*(1+pic_d*total$prec_mm) * total$Urban;  ## Normal
  input$D_pre = Be_devpre*prec_mm_in^pic_d_pre2 * total$`Urban pre-1980`;  ## Normal
  input$D_post = Be_devpost*prec_mm_in^pic_d_post2 * (total$Urban_80_2000+total$`Urban post-2000`);  ##Pre and post
  #D = Be_d*(1+Bp[2]*av_prec) * D.lc.pre + Be_d*(1+Bp[4]*av_prec) * D.lc.post;  ##Pre and post
  input$W = Be_w*prec_mm_in^pic_w2 * total$Undeveloped;
  
  input$Ch=Be_ch*prec_mm_in^pic_ch2 * total$Chickens;
  input$H=Be_ch*prec_mm_in^pic_h2 * total$Hogs;
  input$Cw=Be_cw*prec_mm_in^pic_cw2 * total$Cows;
  input$L=total$Ch+total$H+total$Cw
  
}


############################################
########################################
##############Get expected retention from each subwatershed based on RT and HL
##############  THis had to be extended to make LMS stations go to the lake
####################################################
############################################

#This can be used to adjust source TN in each subwatershed 
############################################################

###Set up matrix to get adjust RT and HL from LMS to lake values for all subwsds82-2017
#c=seq(from=1982,to=2017,by=1)
col.rt=5
rows.rt=length(unique(total[,"Siteno"]))

####Basic matrix to use
subwsd.ret=as.data.frame(matrix(0,nrow=rows.rt,ncol=col.rt))   ####Initialized matrix
subwsd.ret[,1]<-unique(total[,"Siteno"]);    ###LMS names....CORRECT SITEMATCH NAMES
colnames(subwsd.ret)[1:col.rt]<- c("Subwsd","RT_LMS","HL_LMS","RT_lake","HL_lake")

subwsd.ret$RT_LMS=total$stream_rt[match(subwsd.ret$Subwsd,total$Siteno)]
subwsd.ret$HL_LMS=total$res_hydload[match(subwsd.ret$Subwsd,total$Siteno)]
#subwsd.ret$LMS=total$ [match(subwsd.ret$Subwsd,total$Siteno)]

######Import adjustments to 
LMS.lake= read.csv("/My Drive/Post-Doc 2019/Jordan Lake Watershed Model/Longest Paths Subwatersheds/LMS_lake_match_RT_HL.csv",header=TRUE,na.strings = "-9999")
land_match = read.csv("/My Drive/Post-Doc 2019/Jordan Lake Watershed Model/STAN models/Posteriors/Plots/2019-09-27_TN_yearly_3chain_10000/Landcover2019-09-27.csv",header=TRUE,na.strings = "-9999")
land_match <- land_match[land_match[,3] %in% 2008,]  ####This limits land to only 1994-2017 values

names(land_match)
land_match=land_match[,c(2,13:14)]

subwsd.ret$LMS=land_match$new[match(subwsd.ret$Subwsd,land_match$Siteno)]

####Add additional RT and HL for each subwatershed
subwsd.ret$RT_lake=LMS.lake$RT_lake[match(subwsd.ret$LMS,LMS.lake$Siteno)]
subwsd.ret$HL_lake=LMS.lake$HL_lake[match(subwsd.ret$LMS,LMS.lake$Siteno)]

subwsd.ret[,c(2:5)][is.na(subwsd.ret[,c(2:5)])] <- 0 ##Gets rid of NAs in this column for Dischargers
subwsd.ret[,c(5)][(subwsd.ret[,c(5)])==0] <- 999999 ##Gets rid of NAs in this column for Dischargers

######Join RT and HL for each subwatershed
subwsd.ret$RT_total=subwsd.ret$RT_LMS+subwsd.ret$RT_lake
subwsd.ret$HL_total=0
get_weighted_eff_k=function(eff_k){
  #eff_k=sub$mean_eff_k
  Eff_k_denom=prod(eff_k)
  eff_k_num_coeff=Eff_k_denom/eff_k
  Eff_k_num=sum(eff_k_num_coeff)
  Eff_k=Eff_k_denom/(Eff_k_num)
  #Eff_retention=1-exp(-mass_transfer/Eff_k)
  return(Eff_k)
}

#####Loop to join HL for each subwatershed
for (i in 1:length(subwsd.ret[,1])){
  #i=1
  dt=subwsd.ret[i,]
  k=cbind(dt[3],dt[5])
  k2=get_weighted_eff_k(k)
  if (k2>450000){
    k2=999999
  }
  subwsd.ret[i,8]=k2
  }

######Loop to get stream and lake retention expected for each HUC
#########################################################
#Sn=.03
subwsd.ret$retained_stream=0
subwsd.ret$retained_reservoir=0
subwsd.ret$retained_total=0

get_weighted_retention=function(eff_k,mass_transfer){
  Eff_k_denom=prod(eff_k)
  eff_k_num_coeff=Eff_k_denom/eff_k
  Eff_k_num=mass_transfer*sum(eff_k_num_coeff)
  Eff_k=Eff_k_denom/(Eff_k_num/mass_transfer)
  Eff_retention=1-exp(-mass_transfer/Eff_k)
  return(Eff_retention)
}

for (i in 1:length(subwsd.ret[,1])){
  #i=1
  dt=subwsd.ret[i,]
  subwsd.ret[i,9]=exp(-Sn*dt$RT_total)     #####Retained excluding stream losses
  subwsd.ret[i,10]=exp(-Sn2/dt$HL_total)   #####Retained excluding reservoir losses
  subwsd.ret[i,11]=subwsd.ret[i,9]*subwsd.ret[i,10]   #####Retained total
  
}

###########################
#########GEt mean aggregate of TN load by subwatershed - EXEC SUMMARY QUESTION
###########################
total.TN=total
total.TN$TNload=total$A+total$D_pre+total$D_post+total$W+total$L+total$Dischargers
names(total.TN)
sub.TN=total.TN[,c(1,25)] %>%
  group_by(Siteno) %>%
  summarise_each(funs(mean(., na.rm = TRUE)))
sub.TN=as.data.frame(sub.TN)

subwsd.ret$TN.mean=sub.TN$TNload[match(subwsd.ret$Subwsd,sub.TN$Siteno)]
subwsd.ret$TN.mean.lake=subwsd.ret$retained_total*subwsd.ret$TN.mean
subwsd.ret$basin=ungauge_match$Main_Watershed[match(subwsd.ret$Subwsd,ungauge_match$SiteID_NWALT2)]
subwsd.ret.JL=subwsd.ret[subwsd.ret$basin!="FL",]

############For EXECUTIVE SUMMARY########
Mean.ret=sum(subwsd.ret.JL$TN.mean.lake)/sum(subwsd.ret.JL$TN.mean)

t=Sys.Date()
#write.csv(subwsd.ret,paste("Q:/Shared drives/Jordan Lake Modeling/WatershedModeling/Manuscript/Subwatershed Results/HUC_stream retention_rates_TN_prepost_power",t,".csv",sep=""))
summary(subwsd.ret$retained_stream)
########################
#########Clean up total dataset   NOT NEEDED FOR POWER FUNCTION PICS
########################

##GEt rid of negative values 
names(total)
#total[,c(16:23)][(total[,c(16:23)])<0] <- 0 ##Gets rid of NAs in this column for Dischargers
total[,c(16:24)][(total[,c(16:24)])<0] <- 0 ##Gets rid of NAs in this column for Dischargers

##################
###Get total TN with and without stream retention and with/without dischargers
###############

total$TN_total=total$A+total$D_pre+total$D_post+total$W+total$Ch+total$H+total$Cw+total$Discharger
total$TN_total_nodischargers=total$A+total$D_pre+total$D_post+total$W+total$Ch+total$H+total$Cw

####Retention based on each year
total$RT_total=subwsd.ret$RT_total[match(total$Siteno,subwsd.ret$Subwsd)]
total$HL_total=subwsd.ret$HL_total[match(total$Siteno,subwsd.ret$Subwsd)]
###Get mean retention
total$retention_total_mean=subwsd.ret$retained_total[match(total$Siteno,subwsd.ret$Subwsd)]

######Adjust RT and HL based on PIC Q
total$RT_adj=total$RT_total/(1+PIC_q*total$prec_mm)
total$HL_adj=total$HL_total*(1+PIC_q*total$prec_mm)
total$Retention_adj=exp(-Sn*total$RT_adj)* exp(-Sn2/total$HL_adj) ##Flow adjusted retention
#Retention_adj=exp(-Sn*4.72)* exp(-Sn2/10.2) ##Flow adjusted retention

total$TN_total_lake=total$TN_total*total$Retention_adj
total$TN_total_nodischargers_lake=total$TN_total_nodischargers*total$Retention_adj
names(total)
total=total[,c(1:26,33:34,27:32)]

####################################
############Get TN load for predictions for normal, high, low precipitation values
######################################
###############This is done for all years. SHould be reported for just one year
prec.mean=0
prec.10=quantile(total$prec_mm,.1,na.rm=TRUE)
prec.90=quantile(total$prec_mm,.9,na.rm=TRUE)

prec.mean.p=1
prec.10.p=quantile(prec_mm_in,.1,na.rm=TRUE)
prec.90.p=quantile(prec_mm_in,.9,na.rm=TRUE)

###Function to get predictions for different precipitation values

precipitation_TN=function(input,prec){
  #input=total; prec=0   ######Set conditions for prec.
  input$prec_mm=prec
  input$A.p=  (Be_a*(1+pic_a*input$prec_mm) * input$Agriculture);
  input$D.p = Be_d*(1+pic_d*input$prec_mm) * input$Urban;  ## Normal
  #D = Be_d*(1+Bp[2]*av_prec) * D.lc.pre + Be_d*(1+Bp[4]*av_prec) * D.lc.post;  ##Pre and post
  input$W.p = Be_w*(1+pic_w*input$prec_mm) * input$Undeveloped;
  
  input$Ch.p=Be_ch*(1+pic_ch*input$prec_mm) * input$Chickens;
  input$H.p=Be_ch*(1+pic_h*input$prec_mm) * input$Hogs;
  input$Cw.p=Be_cw*(1+pic_cw*input$prec_mm) * input$Cows;
  input$L.p=input$Ch+input$H+input$Cw

  input$TN_total=input$A.p+input$D.p+input$W.p+input$Ch.p+input$H.p+input$Cw.p+input$Discharger
  input$TN_total_nodischargers=input$A.p+input$D.p+input$W.p+input$Ch.p+input$H.p+input$Cw.p

  input$TN_total_lake=input$TN_total*input$retention_total
  input$TN_total_nodischargers_lake=input$TN_total_nodischargers*input$retention_total
  names(input)
  return(input[,c(1:2,24:31)])
   }

precipitation_TN_prepost=function(input,prec){
  #input=total; prec=0   ######Set conditions for prec.
  input$prec_mm=prec
  input$A.p=  (Be_a*(1+pic_a*input$prec_mm) * input$Agriculture);
  #input$D.p = Be_d*(1+pic_d*input$prec_mm) * input$Urban;  ## Normal
  input$D.pre = Be_devpre*(1+pic_d_pre*input$prec_mm) * input$`Urban pre-1980`;  ## Normal
  input$D.post = Be_devpost*(1+pic_d_post*input$prec_mm) * (input$Urban_80_2000+input$`Urban post-2000`);  ## Normal
  #D = Be_d*(1+Bp[2]*av_prec) * D.lc.pre + Be_d*(1+Bp[4]*av_prec) * D.lc.post;  ##Pre and post
  input$W.p = Be_w*(1+pic_w*input$prec_mm) * input$Undeveloped;
  
  input$Ch.p=Be_ch*(1+pic_ch*input$prec_mm) * input$Chickens;
  input$H.p=Be_ch*(1+pic_h*input$prec_mm) * input$Hogs;
  input$Cw.p=Be_cw*(1+pic_cw*input$prec_mm) * input$Cows;
  input$L.p=input$Ch+input$H+input$Cw
  
  input$TN_total=input$A.p+input$D.pre+input$D.post+input$W.p+input$Ch.p+input$H.p+input$Cw.p+input$Discharger
  input$TN_total_nodischargers=input$A.p+input$D.pre+input$D.post+input$W.p+input$Ch.p+input$H.p+input$Cw.p
  
  input$TN_total_lake=input$TN_total*input$retention_total
  input$TN_total_nodischargers_lake=input$TN_total_nodischargers*input$retention_total
  names(input)
  return(input[,c(1:2,24:34)])
}

precipitation_TN_prepost_power_2017=function(input,prec){
  #input=total; prec=1
  input=input[input$Year==2017,]
  input$A.p=  (Be_a*prec^pic_a2 * input$Agriculture);
  #input$D = Be_d*(1+pic_d*total$prec_mm) * inputl$Urban;  ## Normal
  input$D.pre = Be_devpre*prec^pic_d_pre2 * input$`Urban pre-1980`;  ## Normal
  input$D.post = Be_devpost*prec^pic_d_post2 * (input$Urban_80_2000+input$`Urban post-2000`);  ##Pre and post
  #D = Be_d*(1+Bp[2]*av_prec) * D.lc.pre + Be_d*(1+Bp[4]*av_prec) * D.lc.post;  ##Pre and post
  input$W.p = Be_w*prec^pic_w2 * input$Undeveloped;
  
  input$Ch.p=Be_ch*prec^pic_ch2 * input$Chickens;
  input$H.p=Be_ch*prec^pic_h2 * input$Hogs;
  input$Cw.p=Be_cw*prec^pic_cw2 * input$Cows;
  input$L.p=input$Ch.p+input$H.p+input$Cw.p
  
  input$TN_total=input$A.p+input$D.pre+input$D.post+input$W.p+input$L.p+input$Discharger
  input$TN_total_nodischargers=input$A.p+input$D.pre+input$D.post+input$W.p+input$L.p
  
  #####Adj RT and HL for precipitation
  #####Convert prec_mm to prec/mean(prec)
  #mean=1179.6;stdev=171.8
  scale=(prec*1179.6-1179.6)/171.8

  input$RT_adj=input$RT_total/(1+PIC_q*scale)
  input$HL_adj=input$HL_total*(1+PIC_q*scale)
  input$Retention_adj=exp(-Sn*input$RT_adj)* exp(-Sn2/input$HL_adj) ##Flow adjusted retention
  
  input$TN_total_lake=input$TN_total*input$Retention_adj
  input$TN_total_nodischargers_lake=input$TN_total_nodischargers*input$Retention_adj
  names(input)
  return(input)
}


#p.mean=precipitation_TN(total,prec.mean)
#p.10=precipitation_TN(total,prec.10)
#p.90=precipitation_TN(total,prec.90)

p.mean=precipitation_TN_prepost(total,prec.mean)
p.10=precipitation_TN_prepost(total,prec.10)
p.90=precipitation_TN_prepost(total,prec.90)

p.mean.p=precipitation_TN_prepost_power_2017(total,prec.mean.p)
p.10.p=precipitation_TN_prepost_power_2017(total,prec.10.p)
p.90.p=precipitation_TN_prepost_power_2017(total,prec.90.p)


########################################################
######Loop to get time series of TN loading for each HUC
########################################################
years=seq(1994,2017,1)
col=length(years)+4    #####To add normalized precipitation values
rows=length(unique(total$Siteno))
results=as.data.frame(matrix(0,nrow=rows,ncol=col+length(years)))   ####Initialized matrix
results[,1]<-unique(total$Siteno);    ###County names
colnames(results)<- c("Siteno","Prec_10%_2017","Prec_mean_2017","Prec_90%_2017",years,years)

results.no.dsch=results  ###Get results without dischargers
results.ret=results  ###Get results with stream retention
results.no.dsch.ret=results  ###Get results without dischargers with stream retention


########This loop runs depending on if using PIC_Q or PIC_Q
for (i in 1:(col-4)){
  #i=1;
  yr=years[i]
  tot=total[total$Year==yr,]   ###Predictions with real rainfall
  #tot.m=p.mean.p[p.mean.p$Year==yr,] ###Predictions with mean rainfall
  #tot.10=p.10.p[p.10.p$Year==yr,]  ###Predictions with 10% rainfall
  #tot.90=p.90.p[p.90.p$Year==yr,]  ###Predictions with 90% rainfall
  
  for (j in 1:rows){
    #j=1
    site=as.character(results[j,1])  
    tot2=tot[tot$Siteno==site,]
    tot2.m=p.mean.p[p.mean.p$Siteno==site,]
    tot2.10=p.10.p[p.10.p$Siteno==site,]
    tot2.90=p.90.p[p.90.p$Siteno==site,]
    ###Results for actual predictions
    results[j,i+4]=tot2$TN_total
    results.no.dsch[j,i+4]=tot2$TN_total_nodischargers
    results.ret[j,i+4]=tot2$TN_total
    results.no.dsch.ret[j,i+4]=tot2$TN_total_nodischargers_lake
    ####Get retention rates
    results[j,i+4+24]=tot2$Retention_adj
    results.no.dsch[j,i+4+24]=tot2$Retention_adj
    results.ret[j,i+4+24]=tot2$Retention_adj
    results.no.dsch.ret[j,i+4+24]=tot2$Retention_adj###Results for mean prec
    ###Results for mean prec
    results[j,3]=tot2.m$TN_total
    results.no.dsch[j,3]=tot2.m$TN_total_nodischargers
    results.ret[j,3]=tot2.m$TN_total_lake
    results.no.dsch.ret[j,3]=tot2.m$TN_total_nodischargers_lake
    ###Results for 10% prec
    results[j,2]=tot2.10$TN_total
    results.no.dsch[j,2]=tot2.10$TN_total_nodischargers
    results.ret[j,2]=tot2.10$TN_total_lake
    results.no.dsch.ret[j,2]=tot2.10$TN_total_nodischargers_lake
    ###Results for 90% prec
    results[j,4]=tot2.90$TN_total
    results.no.dsch[j,4]=tot2.90$TN_total_nodischargers
    results.ret[j,4]=tot2.90$TN_total_lake
    results.no.dsch.ret[j,4]=tot2.90$TN_total_nodischargers_lake
    
  }
}  ###Old loop for linear....using retention from subwsd.ret


###New loop.. gets retention with PIC_Q

####For loop for retention with PIC_Q2 (square root for Velocity)
# 
#   for (i in 1:(col-4)){
#     #i=14; ###2007
#     yr=years[i]
#     tot=total[total$Year==yr,]   ###Predictions with real rainfall
#     tot.m=p.mean[p.mean$Year==yr,] ###Predictions with mean rainfall
#     tot.10=p.10[p.10$Year==yr,]  ###Predictions with 10% rainfall
#     tot.90=p.90[p.90$Year==yr,]  ###Predictions with 90% rainfall
#     
#     subwsd.ret[i,9]=exp(-Sn*dt$RT_total)     #####Retained excluding stream losses
#     subwsd.ret[i,10]=exp(-Sn2/dt$HL_total)   #####Retained excluding reservoir losses
#     
#     
#     for (j in 1:rows){
#       #j=10
#       site=as.character(results[j,1])  
#       tot2=tot[tot$Siteno==site,]
#       tot2.m=tot.m[tot.m$Siteno==site,]
#       tot2.10=tot.10[tot.10$Siteno==site,]
#       tot2.90=tot.90[tot.90$Siteno==site,]
#       ###Results for actual predictions
#       results[j,i+4]=tot2$TN_total
#       results.no.dsch[j,i+4]=tot2$TN_total_nodischargers
#       PIC_adj=1+PIC_q*tot2$prec_mm
#       PIC_adj2=PIC_adj^.5
#       results.ret[j,i+4]=tot2$TN_total*exp(-Sn*tot2$RT_total/PIC_adj2)*exp((-Sn2/PIC_adj)/tot2$HL_total)
#       results.no.dsch.ret[j,i+4]=tot2$TN_total_nodischargers*exp(-Sn*tot2$RT_total/PIC_adj2)*exp((-Sn2/PIC_adj)/tot2$HL_total)
#       ####Get retention rates
#       Ret=exp(-Sn*tot2$RT_total/PIC_adj2)*exp((-Sn2/PIC_adj)/tot2$HL_total)
#       results[j,i+4+24]=Ret
#       results.no.dsch[j,i+4+24]=Ret
#       results.ret[j,i+4+24]=Ret
#       results.no.dsch.ret[j,i+4+24]=Ret###Results for mean prec
#       results[j,3]=tot2.m$TN_total
#       results.no.dsch[j,3]=tot2.m$TN_total_nodischargers
#       results.ret[j,3]=tot2.m$TN_total*exp(-Sn*tot2.m$RT_total/PIC_adj2)*exp((-Sn2/PIC_adj)/tot2.m$HL_total)
#       results.no.dsch.ret[j,3]=tot2.m$TN_total_nodischargers*exp(-Sn*tot2.m$RT_total/PIC_adj2)*exp((-Sn2/PIC_adj)/tot2.m$HL_total)
#       ###Results for 10% prec
#       results[j,2]=tot2.10$TN_total
#       results.no.dsch[j,2]=tot2.10$TN_total_nodischargers
#       results.ret[j,2]=tot2.10$TN_total*exp(-Sn*tot2.10$RT_total/PIC_adj2)*exp((-Sn2/PIC_adj)/tot2.10$HL_total)
#       results.no.dsch.ret[j,2]=tot2.10$TN_total_nodischargers*exp(-Sn*tot2.10$RT_total/PIC_adj2)*exp((-Sn2/PIC_adj)/tot2.10$HL_total)
#       ###Results for 90% prec
#       results[j,4]=tot2.90$TN_total
#       results.no.dsch[j,4]=tot2.90$TN_total_nodischargers
#       results.ret[j,4]=tot2.90$TN_total*exp(-Sn*tot2.90$RT_total/PIC_adj2)*exp((-Sn2/PIC_adj)/tot2.90$HL_total)
#       results.no.dsch.ret[j,4]=tot2.90$TN_total_nodischargers*exp(-Sn*tot2.90$RT_total/PIC_adj2)*exp((-Sn2/PIC_adj)/tot2.90$HL_total)
#       
#     }
#   }  ###New loop.. gets retention with PIC_Q2 
#   
######Need to   #####Fix dishcargers....
########1) summarize by 3 basins for each year.... to graph like before....
#########2) Dump values for each subwsd for each year into excel sheet for Dan
######## 3) Get 1994,2000, 2010,2018 predictions for subwsds for mean/high/low prec

results$basin=ungauge_match$Main_Watershed[match(results$Siteno,ungauge_match$SiteID_NWALT2)]
results.no.dsch$basin=ungauge_match$Main_Watershed[match(results.no.dsch$Siteno,ungauge_match$SiteID_NWALT2)]
results.ret$basin=ungauge_match$Main_Watershed[match(results.ret$Siteno,ungauge_match$SiteID_NWALT2)]
results.no.dsch.ret$basin=ungauge_match$Main_Watershed[match(results.no.dsch.ret$Siteno,ungauge_match$SiteID_NWALT2)]


#####Export data dump for actual TN load by HUC

t=Sys.Date()
write.csv(results,paste("Q:/Shared drives/Jordan Lake Modeling/WatershedModeling/Manuscript/Subwatershed Results/HUC_pred_TN_prepost",t,".csv",sep=""))
write.csv(results.ret,paste("Q:/Shared drives/Jordan Lake Modeling/WatershedModeling/Manuscript/Subwatershed Results/HUC_pred_TN_with_stream_ret_prepost",t,".csv",sep=""))
write.csv(results.no.dsch,paste("Q:/Shared drives/Jordan Lake Modeling/WatershedModeling/Manuscript/Subwatershed Results/HUC_pred_TN_no_dsch_prepost",t,".csv",sep=""))
write.csv(results.no.dsch.ret,paste("Q:/Shared drives/Jordan Lake Modeling/WatershedModeling/Manuscript/Subwatershed Results/HUC_pred_TN_no_dsch_with_stream_ret_prepost",t,".csv",sep=""))

###########################
######Summaries for Executive summary
###########################

####Bind TN that reaches lake with and without dischargers
#exec=cbind(results.ret[,c(1,3)],results.no.dsch.ret[,c(3,29)])  
exec=cbind(results.ret[,c(1,3)],results.no.dsch.ret[,c(3,53)]) ##More columns now
colnames(exec)[3]<-"Prec_mean2017_nodischargers"

######Aggregate to basins
exec2=exec[,c(2:4)] %>%
  group_by(basin) %>%
  summarise_each(funs(sum(., na.rm = TRUE)))

exec2=as.data.frame(exec2)

###Calculate pct of TN load by discharger and by New Hope
exec2$discharger=exec2$Prec_mean_2017-exec2$Prec_mean2017_nodischargers
exec2$pct_discharger=exec2$discharger/exec2$Prec_mean_2017

######Get Jordan Lake total and NEW HOPE pct contribution
JL_totalTN=exec2[2,2]+exec2[3,2]
NewHope_pct=exec2[3,2]/JL_totalTN  ###FOR EXEC SUMMARY

#####Get Discharger pct to Jordan Lake
Pct_disch_JL=(exec2[2,4]+exec2[3,4])/JL_totalTN  ###FOR EXEC SUMMARY

######25% reduction in pt sources, leads to a XX% reduction of TN to lake
Pct_disch_JL_.25=.75*(exec2[2,4]+exec2[3,4])

Pct_disch_JL_.25=Pct_disch_JL_.25 + (exec2[2,3]+exec2[3,3])
Pct_disch_JL_.25=(JL_totalTN-Pct_disch_JL_.25)/JL_totalTN ###FOR EXEC SUMMARY

#####Get pre urbanization to match post- urbanization
#####Need to rewrite this code for pre/post predictions and then
##### go back to the Total dataset to get these values.
#PIC_q=.5

###Changed, but no adjustment to Retention for mean precipitation
precipitation_TN_pretopost_red_prepost=function(input,prec){
  input=total; prec=0   ######Set conditions for prec.
  input$prec_mm=prec
  Be_devpre=Be_devpost
  input$A.p=  (Be_a*(1+pic_a*input$prec_mm) * input$Agriculture);
  #input$D.p = Be_d*(1+pic_d*input$prec_mm) * input$Urban;  ## Normal
  input$D.pre = Be_devpre*(1+pic_d_pre*input$prec_mm) * input$`Urban pre-1980`;  ## Normal
  input$D.post = Be_devpost*(1+pic_d_post*input$prec_mm) * (input$Urban_80_2000+input$`Urban post-2000`);  ## Normal
  input$W.p = Be_w*(1+pic_w*input$prec_mm) * input$Undeveloped;
  
  input$Ch.p= Be_ch*(1+pic_ch*input$prec_mm) * input$Chickens;
  input$H.p= Be_ch*(1+pic_h*input$prec_mm) * input$Hogs;
  input$Cw.p= Be_cw*(1+pic_cw*input$prec_mm) * input$Cows;
  input$L.p=input$Ch+input$H+input$Cw
  
  input$TN_total=input$A.p+input$D.pre+input$D.post+input$W.p+input$Ch.p+input$H.p+input$Cw.p+input$Discharger
  input$TN_total_nodischargers=input$A.p+input$D.pre+input$D.post+input$W.p+input$Ch.p+input$H.p+input$Cw.p
  
  PIC_adj=1+PIC_q*input$prec_mm
  Ret=exp(-Sn*input$RT_total/PIC_adj)*exp((-Sn2/PIC_adj)/input$HL_total)
  
  input$TN_total_lake=input$TN_total*Ret
  input$TN_total_nodischargers_lake=input$TN_total_nodischargers*Ret
  names(input)
  return(input[,c(1:2,24:29)])
}
precipitation_TN_pretopost_red_prepost_power=function(input,prec){
  #input=total; prec=1   ######Set conditions for prec.
  input=input[input$Year==2017,]
  #input$prec_mm=prec
  Be_devpre=Be_devpost
  input$A.p=  (Be_a*(prec^pic_a2) * input$Agriculture);
  #input$D.p = Be_d*(1+pic_d*input$prec_mm) * input$Urban;  ## Normal
  input$D.pre = Be_devpre*(prec^pic_d_pre2) * input$`Urban pre-1980`;  ## Normal
  input$D.post = Be_devpost*(prec^pic_d_post2) * (input$Urban_80_2000+input$`Urban post-2000`);  ## Normal
  input$W.p = Be_w*(prec^pic_w2) * input$Undeveloped;
  
  input$Ch.p= Be_ch*(prec^pic_ch2) * input$Chickens;
  input$H.p= Be_ch*(prec^pic_h2) * input$Hogs;
  input$Cw.p= Be_cw*(prec^pic_cw2) * input$Cows;
  input$L.p=input$Ch+input$H+input$Cw
  
  input$TN_total=input$A.p+input$D.pre+input$D.post+input$W.p+  input$L.p+input$Discharger
  input$TN_total_nodischargers=input$A.p+input$D.pre+input$D.post+input$W.p+input$L.p
  
  #####Adj RT and HL for precipitation
  #####Convert prec_mm to prec/mean(prec)
  #mean=1179.6;stdev=171.8
  scale=(prec*1179.6-1179.6)/171.8
  
  input$RT_adj=input$RT_total/(1+PIC_q*scale)
  input$HL_adj=input$HL_total*(1+PIC_q*scale)
  input$Retention_adj=exp(-Sn*input$RT_adj)* exp(-Sn2/input$HL_adj) ##Flow adjusted retention
  
  input$TN_total_lake=input$TN_total*input$Retention_adj
  input$TN_total_nodischargers_lake=input$TN_total_nodischargers*input$Retention_adj
  names(input)
  return(input[,c(1:2,24:34)])
}

#p.mean_pretopost_red=precipitation_TN_pretopost_red_prepost(total,prec.mean)

p.mean_pretopost_red=precipitation_TN_pretopost_red_prepost_power(total,1)

dev.t=p.mean.p[p.mean.p$Year==2017,]
dev.red=p.mean_pretopost_red[p.mean_pretopost_red$Year==2017,]

dev.join=cbind(dev.t[,c(1,28)],dev.red[,7])

dev.join$basin=ungauge_match$Main_Watershed[match(dev.join$Siteno,ungauge_match$SiteID_NWALT2)]

######Aggregate by basins
dev.join2=dev.join[,c(2:4)] %>%
  group_by(basin) %>%
  summarise_each(funs(sum(., na.rm = TRUE)))

JL_total_dev_reduction=(dev.join2[2,2]-dev.join2[2,3]+(dev.join2[3,2]-dev.join2[3,3]))
###FOR EXEC SUMMARY
Dev_reduction_pretopost=JL_total_dev_reduction/JL_totalTN ###FOR EXEC SUMMARY

############################################################
##### Get livestock and ag to be reduced 25%
###THese values are slightly off from the ones above??? CHeck
precipitation_TN_ag_red=function(input,prec,reduction){
  #input=total; prec=0   ######Set conditions for prec.
  input$prec_mm=prec
  input$A.p=  (1-reduction)*(Be_a*(1+pic_a*input$prec_mm) * input$Agriculture);
  input$D.p = Be_d*(1+pic_d*input$prec_mm) * input$Urban;  ## Normal
  #D = Be_d*(1+Bp[2]*av_prec) * D.lc.pre + Be_d*(1+Bp[4]*av_prec) * D.lc.post;  ##Pre and post
  input$W.p = Be_w*(1+pic_w*input$prec_mm) * input$Undeveloped;
  
  input$Ch.p= (1-reduction)*Be_ch*(1+pic_ch*input$prec_mm) * input$Chickens;
  input$H.p= (1-reduction)*Be_ch*(1+pic_h*input$prec_mm) * input$Hogs;
  input$Cw.p= (1-reduction)*Be_cw*(1+pic_cw*input$prec_mm) * input$Cows;
  input$L.p=input$Ch+input$H+input$Cw
  
  input$TN_total=input$A.p+input$D.p+input$W.p+input$Ch.p+input$H.p+input$Cw.p+input$Discharger
  input$TN_total_nodischargers=input$A.p+input$D.p+input$W.p+input$Ch.p+input$H.p+input$Cw.p
  
  input$TN_total_lake=input$TN_total*input$retention_total
  input$TN_total_nodischargers_lake=input$TN_total_nodischargers*input$retention_total
  names(input)
  return(input[,c(1:2,24:28)])
}
precipitation_TN_ag_red_prepost=function(input,prec,reduction){
  input=total; prec=0   ######Set conditions for prec.
  input$prec_mm=prec
  input$A.p=  (1-reduction)*(Be_a*(1+pic_a*input$prec_mm) * input$Agriculture);
  #input$D.p = Be_d*(1+pic_d*input$prec_mm) * input$Urban;  ## Normal
  input$D.pre = Be_devpre*(1+pic_d_pre*input$prec_mm) * input$`Urban pre-1980`;  ## Normal
  input$D.post = Be_devpost*(1+pic_d_post*input$prec_mm) * (input$Urban_80_2000+input$`Urban post-2000`);  ## Normal
  input$W.p = Be_w*(1+pic_w*input$prec_mm) * input$Undeveloped;
  
  input$Ch.p= (1-reduction)*Be_ch*(1+pic_ch*input$prec_mm) * input$Chickens;
  input$H.p= (1-reduction)*Be_ch*(1+pic_h*input$prec_mm) * input$Hogs;
  input$Cw.p= (1-reduction)*Be_cw*(1+pic_cw*input$prec_mm) * input$Cows;
  input$L.p=input$Ch+input$H+input$Cw
  
  input$TN_total=input$A.p+input$D.pre+input$D.post+input$W.p+input$Ch.p+input$H.p+input$Cw.p+input$Discharger
  input$TN_total_nodischargers=input$A.p+input$D.pre+input$D.post+input$W.p+input$Ch.p+input$H.p+input$Cw.p
  
  input$TN_total_lake=input$TN_total*input$retention_total
  input$TN_total_nodischargers_lake=input$TN_total_nodischargers*input$retention_total
  names(input)
  return(input[,c(1:2,24:29)])
}


precipitation_TN_ag_red_prepost_power=function(input,prec,reduction){
  input=total; prec=1;reduction=.25   ######Set conditions for prec.
  input=input[input$Year==2017,]
  input$A.p=  (1-reduction)*(Be_a*(prec^pic_a2) * input$Agriculture);
  #input$D.p = Be_d*(1+pic_d*input$prec_mm) * input$Urban;  ## Normal
  input$D.pre = Be_devpre*(prec^pic_d_pre2) * input$`Urban pre-1980`;  ## Normal
  input$D.post = Be_devpost*(prec^pic_d_post2) * (input$Urban_80_2000+input$`Urban post-2000`);  ## Normal
  input$W.p = Be_w*(prec^pic_w2) * input$Undeveloped;
  
  input$Ch.p= Be_ch*(prec^pic_ch2) * input$Chickens;
  input$H.p= Be_ch*(prec^pic_h2) * input$Hogs;
  input$Cw.p= Be_cw*(prec^pic_cw2) * input$Cows;
  input$L.p=(1-reduction)*(input$Ch+input$H+input$Cw)
  
  input$TN_total=input$A.p+input$D.pre+input$D.post+input$W.p+input$L.p+input$Discharger
  input$TN_total_nodischargers=input$A.p+input$D.pre+input$D.post+input$W.p+input$L.p
  
  #####Adj RT and HL for precipitation
  #####Convert prec_mm to prec/mean(prec)
  #mean=1179.6;stdev=171.8
  scale=(prec*1179.6-1179.6)/171.8
  
  input$RT_adj=input$RT_total/(1+PIC_q*scale)
  input$HL_adj=input$HL_total*(1+PIC_q*scale)
  input$Retention_adj=exp(-Sn*input$RT_adj)* exp(-Sn2/input$HL_adj) ##Flow adjusted retention
  
  input$TN_total_lake=input$TN_total*input$Retention_adj
  input$TN_total_nodischargers_lake=input$TN_total_nodischargers*input$Retention_adj
  names(input)
  return(input[,c(1:2,24:34)])
}


#precipitation_TN_ag_red_prepost_power(total,1,.25)
p.mean_ag_red=precipitation_TN_ag_red_prepost_power(total,1,.25)

ag.t=p.mean.p[p.mean.p$Year==2017,]
ag.red=p.mean_ag_red[p.mean_ag_red$Year==2017,]
names(ag.t)
ag.join=cbind(ag.t[,c(1,28)],ag.red[,7])

ag.join$basin=ungauge_match$Main_Watershed[match(ag.join$Siteno,ungauge_match$SiteID_NWALT2)]

######Aggregate by basins
ag.join2=ag.join[,c(2:4)] %>%
  group_by(basin) %>%
  summarise_each(funs(sum(., na.rm = TRUE)))

JL_total_ag_reduction=(ag.join2[2,2]-ag.join2[2,3]+(ag.join2[3,2]-ag.join2[3,3]))

Ag_reduction_.25=JL_total_ag_reduction/JL_totalTN ###FOR EXEC SUMMARY

#####5% reduction for a 25% reduction in Ag lands

########################################################
#####Combination scenario...25% dischargers/ag/livestock; pre to post developed
###########################################################

precipitation_TN_combo_red_prepost=function(input,prec,reduction){
  #input=total; prec=0   ######Set conditions for prec.
  input$prec_mm=prec
  Be_devpre=Be_devpost
  input$A.p=  (1-reduction)*(Be_a*(1+pic_a*input$prec_mm) * input$Agriculture);
  #input$D.p = Be_d*(1+pic_d*input$prec_mm) * input$Urban;  ## Normal
  input$D.pre = Be_devpre*(1+pic_d_pre*input$prec_mm) * input$`Urban pre-1980`;  ## Normal
  input$D.post = Be_devpost*(1+pic_d_post*input$prec_mm) * (input$Urban_80_2000+input$`Urban post-2000`);  ## Normal
  input$W.p = Be_w*(1+pic_w*input$prec_mm) * input$Undeveloped;
  
  input$Ch.p= (1-reduction)*Be_ch*(1+pic_ch*input$prec_mm) * input$Chickens;
  input$H.p= (1-reduction)*Be_ch*(1+pic_h*input$prec_mm) * input$Hogs;
  input$Cw.p= (1-reduction)*Be_cw*(1+pic_cw*input$prec_mm) * input$Cows;
  input$L.p=input$Ch+input$H+input$Cw
  
  input$TN_total=input$A.p+input$D.pre+input$D.post+input$W.p+input$Ch.p+input$H.p+input$Cw.p+(1-reduction)*input$Discharger
  input$TN_total_nodischargers=input$A.p+input$D.pre+input$D.post+input$W.p+input$Ch.p+input$H.p+input$Cw.p
  
  input$TN_total_lake=input$TN_total*input$retention_total
  input$TN_total_nodischargers_lake=input$TN_total_nodischargers*input$retention_total
  names(input)
  return(input[,c(1:2,24:29)])
}
precipitation_TN_combo_red_prepost_power=function(input,prec,reduction){
  input=total; prec=1;reduction=.25   ######Set conditions for prec.
  input=input[input$Year==2017,]
  Be_devpre=Be_devpost
  input$A.p=  (1-reduction)*(Be_a*(prec^pic_a2) * input$Agriculture);
  #input$D.p = Be_d*(1+pic_d*input$prec_mm) * input$Urban;  ## Normal
  input$D.pre = Be_devpre*(prec^pic_d_pre2) * input$`Urban pre-1980`;  ## Normal
  input$D.post = Be_devpost*(prec^pic_d_post2) * (input$Urban_80_2000+input$`Urban post-2000`);  ## Normal
  input$W.p = Be_w*(prec^pic_w2) * input$Undeveloped;
  
  input$Ch.p= Be_ch*(prec^pic_ch2) * input$Chickens;
  input$H.p= Be_ch*(prec^pic_h2) * input$Hogs;
  input$Cw.p= Be_cw*(prec^pic_cw2) * input$Cows;
  input$L.p=(1-reduction)*(input$Ch+input$H+input$Cw)
  
  input$TN_total=input$A.p+input$D.pre+input$D.post+input$W.p+input$L.p+ (1-reduction)*input$Discharger
  input$TN_total_nodischargers=input$A.p+input$D.pre+input$D.post+input$W.p+input$L.p
  
  #####Adj RT and HL for precipitation
  #####Convert prec_mm to prec/mean(prec)
  #mean=1179.6;stdev=171.8
  scale=(prec*1179.6-1179.6)/171.8
  
  input$RT_adj=input$RT_total/(1+PIC_q*scale)
  input$HL_adj=input$HL_total*(1+PIC_q*scale)
  input$Retention_adj=exp(-Sn*input$RT_adj)* exp(-Sn2/input$HL_adj) ##Flow adjusted retention
  
  input$TN_total_lake=input$TN_total*input$Retention_adj
  input$TN_total_nodischargers_lake=input$TN_total_nodischargers*input$Retention_adj
  names(input)
  return(input[,c(1:2,24:34)])
}

p.mean_combo_red=precipitation_TN_combo_red_prepost_power(total,prec.mean,.25)

combo.t=p.mean.p[p.mean.p$Year==2017,]
combo.red=p.mean_combo_red[p.mean_combo_red$Year==2017,]
#names(combo.t)
combo.join=cbind(combo.t[,c(1,27)],combo.red[,6])

combo.join$basin=ungauge_match$Main_Watershed[match(combo.join$Siteno,ungauge_match$SiteID_NWALT2)]

######Aggregate by basins
combo.join2=combo.join[,c(2:4)] %>%
  group_by(basin) %>%
  summarise_each(funs(sum(., na.rm = TRUE)))

JL_total_combo=(combo.join2[2,2]-combo.join2[2,3]+(combo.join2[3,2]-combo.join2[3,3]))

###FOR EXEC SUMMARY 
Combo_reduction_.25=JL_total_combo/JL_totalTN ###FOR EXEC SUMMARY

###Combo reduction 30%

############################
#####################Comparing top 25% of years and bottom 25% of years
###########
names(results.ret)
TNtotal.low=results.ret[,c(1,5,8,11,12,16,18,29)]

######Aggregate by basins
TNtotal.low2=TNtotal.low[,c(2:8)] %>%
  group_by(basin) %>%
  summarise_each(funs(sum(., na.rm = TRUE)))

TNtotal.low2=TNtotal.low2[c(2:3),]
TN.low=mean(apply(TNtotal.low2[,c(2:7)],2,sum))

#####High flow years
TNtotal.high=results.ret[,c(1,6,7,9,14,20,24,26,29)]

######Aggregate by basins
TNtotal.high2=TNtotal.high[,c(2:9)] %>%
  group_by(basin) %>%
  summarise_each(funs(sum(., na.rm = TRUE)))

TNtotal.high2=TNtotal.high2[c(2:3),]
TN.high=mean(apply(TNtotal.high2[,c(2:7)],2,sum))

####Ratio of low to high flow years
TNratio=TN.high/TN.low  ###FOR EXEC SUMMARY

##############################
######Basin summaries
##############################
names(results.ret)
names(total)

#####Get only the columns of interest from total
bas=total[,c(15,2,16:19,23:24)]   ####

######NEED to FIX THIS.... TN that reaches the lake!!!!!
names(total)
#bas.ret=total[,c(15,2,25:27,28:29)]  ####
bas.ret=total[,c(15,2,25:28)]  ####

bas[,c(3:8)][(bas[,c(3:8)])<0] <- 0 ##Gets rid of NAs in this column for Dischargers
bas.ret[,c(3:7)][(bas.ret[,c(3:7)])<0] <- 0 ##Gets rid of NAs in this column for Dischargers


####Join all HUCS in the same basin 
bas2=bas %>%
  group_by(basin,Year) %>%
  summarise_each(funs(sum(., na.rm = TRUE)))
bas2=as.data.frame(bas2)


bas.ret2=bas.ret %>%
  group_by(basin,Year) %>%
  summarise_each(funs(sum(., na.rm = TRUE)))
bas.ret2=as.data.frame(bas.ret2)

#####Change column names for ggplot2 
colnames(bas2)[3:8]<-c("Agriculture","Urban, pre-1980","Urban, post-1980","Undeveloped","Livestock","Discharger")
#colnames(bas.ret2)[3:8]<-c("Agriculture","Urban","Undeveloped","Livestock","Discharger")

#####Add the correct discharge data (aggregated above in code)
#bas2.d=merge(bas2,disch2,by=c("basin","Year"),all.x=T,all.y=F); #Seems redundant
#bas2.d=bas2.d[,c(1:6,8)];colnames(bas2.d)[7]="Discharger"

#####Add the correct discharge data with STream retention????? (aggregated above in code)
#basret2.d=merge(bas2,disch2,by=c("basin","Year"),all.x=T,all.y=F);
#basret2.d=bas2.d[,c(1:6,8)];colnames(bas2.d)[7]="Discharger"

#####Melt down the basin data for ggplot2
bas3=melt(bas2,id.vars = c("basin","Year"));
bas.ret3=melt(bas.ret2,id.vars = c("basin","Year"));
#colnames(bas2)[2]<-"Year"; bas2$Year=as.numeric(as.character(bas2$Year))

ttt=bas2$Discharger-bas2.d$Discharger
summary(ttt)


##################Get percents by basin by year as a chart.
###########################
###########Get total for Jordan Lake to do percents
trial.JL=bas3[bas3$basin!="Falls Lake",]
trial.JL=trial.JL[trial.JL$Year> 1993,]
trial.JL2=dcast(trial.JL, variable ~ Year, value.var = "value",fun.aggregate=sum)
means.JL=apply(trial.JL2[,c(2:26)],2,sum)



######################## HAw River Arm
#######################################
trial.HR=bas3[bas3$basin=="Haw River, Jordan Lake",]
trial.HR=trial.HR[trial.HR$Year> 1993,]
trial.HR2=dcast(trial.HR, variable ~ Year, value.var = "value")
means=apply(trial.HR2[,c(2:26)],2,sum)

trial.HR2.pct=trial.HR2
for (i in 2:26){
#i=2
#trial.HR2.pct[,i]=trial.HR2[,i]/sum(trial.HR2[,i]) ###Percent of HAw
trial.HR2.pct[,i]=trial.HR2[,i]/means.JL[(i-1)] ###Percent of Jordan Lake
}
total.pcts=apply(trial.HR2.pct[,c(2:26)],2,sum)
mean(total.pcts[c(9,11,13,15,17,18,21,23)]) ###Normal flow pct of load
mean(total.pcts[c(1,4,7,8,12,14,19,24)])##Low flow pct of load
mean(total.pcts[c(2,3,5,6,10,16,20,22)]) ###High flow pct of load

###Normal years
names(trial.HR2.pct)
trial.HR2.pct.n=trial.HR2.pct[,c(10,12,14,16,18,19,22,24)] ###Normal flow years
norm.means.HR=apply(trial.HR2.pct.n,1,mean)

###Low flow years
names(trial.HR2.pct)
trial.HR2.pct.l=trial.HR2.pct[,c(2,5,8,9,13,15,20,25)] ###Low flow years
low.means.HR=apply(trial.HR2.pct.l,1,mean)

###High flow years
names(trial.HR2.pct)
trial.HR2.pct.h=trial.HR2.pct[,c(3,4,6,7,11,17,21,23)] ###High flow years
high.means.HR=apply(trial.HR2.pct.h,1,mean)

######################## New Hope Arm
#######################################
trial.NH=bas3[bas3$basin=="New Hope Creek, Jordan Lake",]
trial.NH=trial.NH[trial.NH$Year> 1993,]
trial.NH2=dcast(trial.NH, variable ~ Year, value.var = "value")
means=apply(trial.NH2[,c(2:26)],2,sum)

trial.NH2.pct=trial.NH2
for (i in 2:26){
  #i=2
  #trial.NH2.pct[,i]=trial.NH2[,i]/sum(trial.NH2[,i])###Percent of NEw Hope
  trial.NH2.pct[,i]=trial.NH2[,i]/means.JL[(i-1)] ###Percent of Jordan Lake
  }
total.pcts=apply(trial.NH2.pct[,c(2:26)],2,sum)

###Normal years
names(trial.NH2.pct)
trial.NH2.pct.n=trial.NH2.pct[,c(10,12,14,16,18,19,22,24)] ###Normal flow years
norm.means.NH=apply(trial.NH2.pct.n,1,mean)

###Low flow years
names(trial.NH2.pct)
trial.NH2.pct.l=trial.NH2.pct[,c(2,5,8,9,13,15,20,25)] ###Low flow years
low.means.NH=apply(trial.NH2.pct.l,1,mean)

###High flow years
names(trial.NH2.pct)
trial.NH2.pct.h=trial.NH2.pct[,c(3,4,6,7,11,17,21,23)] ###High flow years
high.means.NH=apply(trial.NH2.pct.h,1,mean)

#######################################
trial.FL=bas3[bas3$basin=="Falls Lake",]
trial.FL=trial.FL[trial.FL$Year> 1993,]
trial.FL2=dcast(trial.FL, variable ~ Year, value.var = "value")
means=apply(trial.FL2[,c(2:26)],2,sum)

trial.FL2.pct=trial.FL2
for (i in 2:26){
  #i=2
  trial.FL2.pct[,i]=trial.FL2[,i]/sum(trial.FL2[,i])
}
total.pcts=apply(trial.FL2.pct[,c(2:26)],2,sum)

###Normal years
names(trial.FL2.pct)
trial.FL2.pct.n=trial.FL2.pct[,c(10,12,14,16,18,19,22,24)] ###Normal flow years
norm.means.FL=apply(trial.FL2.pct.n,1,mean)

###Low flow years
names(trial.FL2.pct)
trial.FL2.pct.l=trial.FL2.pct[,c(2,5,8,9,13,15,20,25)] ###Low flow years
low.means.FL=apply(trial.FL2.pct.l,1,mean)

###High flow years
names(trial.FL2.pct)
trial.FL2.pct.h=trial.FL2.pct[,c(3,4,6,7,11,17,21,23)] ###High flow years
high.means.FL=apply(trial.FL2.pct.h,1,mean)

##################################
#######Join all values together
#################################
col=10    
rows=length(trial.HR2.pct[,1])
source.pct=as.data.frame(matrix(0,nrow=rows,ncol=col))   ####Initialized matrix
source.pct[,1]<-trial.HR2.pct[,1];    ###Source names
colnames(source.pct)<- c("Source","HR_low","HR_normal","HR_high","NH_low","NH_normal","NH_high","FL_low","FL_normal","FL_high")
source.pct[,2]=low.means.HR
source.pct[,3]=norm.means.HR
source.pct[,4]=high.means.HR
source.pct[,5]=low.means.NH
source.pct[,6]=norm.means.NH
source.pct[,7]=high.means.NH
source.pct[,8]=low.means.FL
source.pct[,9]=norm.means.FL
source.pct[,10]=high.means.FL

t=Sys.Date()
write.csv(source.pct,paste("Q:/Shared drives/Jordan Lake Modeling/WatershedModeling/Manuscript/Subwatershed Results/BAR CHART_as table_TN",t,".csv",sep=""))

##################Get ration of TN for high and low flow years
trial.f=bas3[bas3$basin!="Falls Lake",]
trial.f=trial.f[,c(2:4)]
trial.f=trial.f[trial.f$Year> 1993,]
trial.f2=dcast(trial.f, variable ~ Year, value.var = "value",fun.aggregate = sum)

trial.f.low=trial.f2[,c(2,5,8,9,13,15,20,25)] ###Low flow years
trial.f.norm=trial.f2[,c(10,12,14,16,18,19,22,24)] ###Normal flow years
trial.f.high=trial.f2[,c(3,4,6,7,11,17,21,23)] ###High flow years


means.low=mean(apply(trial.f.low[,c(2:8)],2,sum))
means.norm=mean(apply(trial.f.norm[,c(2:8)],2,sum))
means.high=mean(apply(trial.f.high[,c(2:8)],2,sum))

#####For executive summary
ratio=means.high/means.low

###Normal years
names(trial.HR2.pct)
trial.HR2.pct.n=trial.HR2.pct[,c(10,12,14,16,18,19,22,24)] ###Normal flow years
norm.means.HR=apply(trial.HR2.pct.n,1,mean)


####################################
######Graphing functions
#####################################
library(ggplot2)

####Get precipitation data to add to secondary axis
####
prec = read.csv("/Shared drives/Jordan Lake Modeling/WatershedModeling/Data/PRISM/Processed Data/Yearly_Prec2019-10-17.csv",header=TRUE,na.strings = "-9999")
pct66=quantile(prec$Prec,probs=c(.10,.5,.90))
prec$basin=ungauge_match$Main_Watershed[match(prec$Siteno,ungauge_match$SiteID_NWALT2)]

##Get columns needed and aggregate to the basin... get mean precipitation by year for the basin
prec=prec[,c(3:5)]

prec.2=prec %>%
  group_by(basin,Year) %>%
  summarise_each(funs(mean(., na.rm = TRUE)))
prec.2=as.data.frame(prec.2)

######Rename basins in prec dataset
levels(prec.2$basin)[levels(prec.2$basin)=="HR"] <- "Haw River, Jordan Lake"
levels(prec.2$basin)[levels(prec.2$basin)=="NHC"] <- "New Hope Creek, Jordan Lake"
levels(prec.2$basin)[levels(prec.2$basin)=="FL"] <- "Falls Lake"

####Merge Mean precipitation to the dataset for graphing 
bar.pr=merge(bas3,prec.2,by=c("basin","Year"),all.x=T,all.y=F)
bar.pr=bar.pr[complete.cases(bar.pr), ]

#####################Getting year with 10%, 50%, and 90% precipitation by basin
bas.total=bas3[,c(1,2,4)] %>%
  group_by(basin,Year) %>%
  summarise_each(funs(sum(., na.rm = TRUE)))

dummy2 <- data.frame(X = c("Haw River Arm, JL", "Haw River Arm, JL","Haw River Arm, JL"
                           ), Z = c(1206581,1702608,2562808))



###Exporting of Bar graph time series
t=Sys.Date()
png(filename=paste("Q:/Shared drives/Jordan Lake Modeling/WatershedModeling/Manuscript/Figures/Bar charts/TN 3basins_stream_ret_prepost",t,".png",sep=""), 
    units="cm",width=20,height=16,res=300)

theme_bw()
bar=bar.pr %>% filter(Year > 1994 & Year < 2018) %>% 
  ggplot(., aes(x = Year, y = value / 100000, fill = factor(variable,levels=c("Livestock","Undeveloped","Agriculture","Urban, pre-1980","Urban, post-1980","Discharger")))) + 
  geom_bar(stat = "identity")+
  facet_wrap(~basin,scales="free",ncol=1)+
  
  ##Precipitation line
  #geom_line(data=bar.pr, aes(x=Year, y=Prec/100 ), color="blue")+
  #geom_hline(data = dummy2, aes(yintercept = Z))+
  #scale_y_continuous(sec.axis = sec_axis(~./max(bar.pr$Prec/100)))+
  theme_bw()+
  theme(panel.grid.major.x = element_blank(),panel.grid.minor.x = element_blank())+
  #theme_classic()+
  scale_fill_manual(values = c("#5e3c99","#b2abd2","khaki2","#fdb863","#e66101","#9ebcda"))+
  theme(legend.title = element_blank())+    #scale_x_continuous(breaks=seq(0,82,1))+
  #ggtitle("TN sources by basin (no retention)")+
  #ylab("TN x 10 (kg/yr)")+xlab("")+
  ylab(bquote('TN x '* ~ 10^5*' (kg/yr)'))+xlab("")+
  theme(legend.title = element_blank(),axis.text.x = element_text(size=10),)+
  scale_x_continuous(breaks=c(1995,2000,2005,2010,2015))
  
bar
dev.off()

#####Bar chart without Falls Lake

png(filename=paste("Q:/Shared drives/Jordan Lake Modeling/WatershedModeling/Manuscript/Figures/Bar charts/TN 3basins_stream_ret_prepost_JLonly",t,".png",sep=""), 
    units="cm",width=20,height=16,res=300)

bar.pr.JL=bar.pr[bar.pr$basin!="Falls Lake",]

theme_bw()
bar=bar.pr.JL %>% filter(Year > 1994 & Year < 2018) %>% 
  ggplot(., aes(x = Year, y = value / 100000, fill = factor(variable,levels=c("Livestock","Undeveloped","Agriculture","Urban, pre-1980","Urban, post-1980","Discharger")))) + 
  geom_bar(stat = "identity")+
  
  ##Precipitation line
  #geom_line(data=bar.pr, aes(x=Year, y=Prec/100 ), color="blue")+
  #geom_hline(data = dummy2, aes(yintercept = Z))+
  #scale_y_continuous(sec.axis = sec_axis(~./max(bar.pr$Prec/100)))+
  theme_bw(base_size = 14)+
  theme(panel.grid.major.x = element_blank(),panel.grid.minor.x = element_blank())+
  #theme_classic()+
  facet_wrap(~basin,scales="free",ncol=1)+
  scale_fill_manual(values = c("#5e3c99","#b2abd2","khaki2","#fdb863","#e66101","#9ebcda"))+
  theme(legend.title = element_blank())+    #scale_x_continuous(breaks=seq(0,82,1))+
  #ggtitle("TN sources by basin (no retention)")+
  #ylab("TN x 10 (kg/yr)")+xlab("")+
  ylab(bquote('TN x '* ~ 10^5*' (kg/yr)'))+xlab("")+
  #theme(legend.title = element_blank(),element_text(size=14),axis.text.x = element_text(size=16),strip.text.x = element_text(size = 12))+
  theme(legend.title = element_blank())+
  scale_x_continuous(breaks=c(1995,2000,2005,2010,2015))

bar
dev.off()

#############################################################
################With error bars from random effect
#############################################################

#####Melt down the basin data for ggplot2
Total=bas2$Agriculture+bas2$`Urban, pre-1980`+bas2$`Urban, post-1980`+bas2$Undeveloped+bas2$Livestock+bas2$Discharger
bas3.e=melt(bas2,id.vars = c("basin","Year"));

####Merge Mean precipitation to the dataset for graphing 
bas3$Total=
bar.pr.e=merge(bas3,prec.2,by=c("basin","Year"),all.x=T,all.y=F)
bar.pr.e=bar.pr.e[complete.cases(bar.pr), ]


###Exporting of Bar graph time series
t=Sys.Date()
png(filename=paste("Q:/Shared drives/Jordan Lake Modeling/WatershedModeling/Manuscript/Figures/Bar charts/TN 3basins_stream_ret_prepost_errorbars",t,".png",sep=""), 
    units="cm",width=20,height=16,res=300)

theme_bw()
bar=bar.pr %>% filter(Year > 1994 & Year < 2018) %>% 
  ggplot(., aes(x = Year, y = value / 100000, fill = factor(variable,levels=c("Livestock","Undeveloped","Agriculture","Urban, pre-1980","Urban, post-1980","Discharger")))) + 
  geom_bar(stat = "identity")+
  
  ##Precipitation line
  #geom_line(data=bar.pr, aes(x=Year, y=Prec/100 ), color="blue")+
  #geom_hline(data = dummy2, aes(yintercept = Z))+
  #scale_y_continuous(sec.axis = sec_axis(~./max(bar.pr$Prec/100)))+
  theme_bw()+
  theme(panel.grid.major.x = element_blank(),panel.grid.minor.x = element_blank())+
  #theme_classic()+
  facet_wrap(~basin,scales="free",ncol=1)+
  scale_fill_manual(values = c("#5e3c99","#b2abd2","khaki2","#fdb863","#e66101","#5e3c99"))+
  theme(legend.title = element_blank())+    #scale_x_continuous(breaks=seq(0,82,1))+
  #ggtitle("TN sources by basin (no retention)")+
  #ylab("TN x 10 (kg/yr)")+xlab("")+
  ylab(bquote('TN x '* ~ 10^5*' (kg/yr)'))+xlab("")+
  theme(legend.title = element_blank(),axis.text.x = element_text(size=10),)+
  scale_x_continuous(breaks=c(1995,2000,2005,2010,2015))

bar
dev.off()


######2nd version 

t=ggplot(bar.pr)  + 
  geom_bar(aes(x=Year, y=value / 100000),stat="identity", fill = factor(variable,levels=c("Livestock","Undeveloped","Agriculture","Urban","Discharger")
  #geom_line(aes(x=Year, y=Rate*max(df$Response)),stat="identity")+
t
  geom_bar(stat = "identity")+
  +geom_line(data=bar.pr, aes(x=Year, y=Prec/50 ), color="blue")+
  theme_bw()+
  theme(panel.grid.major.x = element_blank(),panel.grid.minor.x = element_blank())+
  #theme_classic()+
  facet_wrap(~basin,scales="free",ncol=1)+
  scale_fill_manual(values = c("#5e3c99","#b2abd2","khaki2","#fdb863","#e66101"))+
  theme(legend.title = element_blank())+    #scale_x_continuous(breaks=seq(0,82,1))+
  ggtitle("TN sources by basin (no retention)")+
  #ylab("TN x 10 (kg/yr)")+xlab("")+
  ylab(bquote('TN x '* ~ 10^5*' (kg/yr)'))+xlab("")+
  theme(legend.title = element_blank(),axis.text.x = element_text(size=10),)+
  scale_x_continuous(breaks=c(1995,2000,2005,2010,2015))
bar
bar+geom_line(data=prec.2, aes(x=Year, y=Prec/50 ), color="blue")+
  scale_y_continuous(limits=c(0, 1200),
 sec.axis = sec_axis(~ . *scale_of_the_new_axis, name = "name of the new axis"))


p <- p + geom_line(aes(y = rel_hum/5, colour = "Humidity"))

# now adding the secondary axis, following the example in the help file ?scale_y_continuous
# and, very important, reverting the above transformation
p <- p + scale_y_continuous(sec.axis = sec_axis(~.*5, name = "Relative humidity [%]"))


##########
#####Get bar chart source vs. stream retention
##########

TN.1995<- bas3[bas3[,2] %in% 1995,]
TN.1995.ret<- bas.ret3[bas.ret3[,2] %in% 1995,]

TN.2005<- bas3[bas3[,2] %in% 2005,]
TN.2005.ret<- bas.ret3[bas.ret3[,2] %in% 2005,]

TN.2015<- bas3[bas3[,2] %in% 2015,]
TN.2015.ret<- bas.ret3[bas.ret3[,2] %in% 2015,]

######Join years of interest
TN.yr=rbind(TN.1995,TN.2005,TN.2015)
TN.yr.ret=rbind(TN.1995.ret,TN.2005.ret,TN.2015.ret)
TN.yr$type="No retention"
TN.yr.ret$type="Stream/Lake retention"

TN.total=rbind(TN.yr,TN.yr.ret)

graph1994=function(data,xlab){
  ggplot(data, aes(x = basin, y = value/ 100000, fill = factor(variable,levels=c("UpstreamTNload","Livestock","Forest","Agriculture","Urban","Discharger")))) + 
    geom_bar(colour="black",stat = "identity")+
    theme_classic()+
    scale_fill_manual(values = c("khaki2", "#e66101","#fdb863","#b2abd2","#5e3c99"))+
    theme(legend.title = element_blank())+
    #scale_fill_viridis(discrete=TRUE,option = "")+
    #scale_x_continuous(breaks=seq(0,82,1))+
    ggtitle("1994 TN sources")+
    ylab(bquote('TN x '* ~ 10^5*' (kg/yr)'))+xlab(paste(xlab))+
    theme(legend.title = element_blank(),axis.text.x = element_text(size=10),) 
  
}
graph3years=function(data,xlab){
  ggplot(data, aes(x = basin, y = value/ 100000, fill = factor(variable,levels=c("Livestock","Undeveloped","Agriculture","Urban","Discharger")))) + 
    geom_bar(stat = "identity")+
    theme_classic()+
    scale_fill_manual(values = c("khaki2", "#e66101","#fdb863","#b2abd2","#5e3c99"))+
    facet_wrap(~Year,scales="free",ncol=1)+
    theme(legend.title = element_blank())+
    #scale_fill_viridis(discrete=TRUE,option = "")+
    #scale_x_continuous(breaks=seq(0,82,1))+
    ggtitle("TN sources")+
    ylab(bquote('TN x '* ~ 10^5*' (kg/yr)'))+xlab(paste(xlab))+
    theme(legend.title = element_blank(),axis.text.x = element_text(size=10),) 
  
}
graph3years_retvsnone=function(data,xlab){
  ggplot(data, aes(x = basin, y = value/ 100000, fill = factor(variable,levels=c("Livestock","Undeveloped","Agriculture","Urban","Discharger")))) + 
    geom_bar(colour="black",stat = "identity")+
    theme_bw()+
    scale_fill_manual(values = c("khaki2", "#e66101","#fdb863","#b2abd2","#5e3c99"))+
    theme(panel.grid.major.x = element_blank(),panel.grid.minor.x = element_blank())+
    #facet_wrap(~Year,scales="free",ncol=2)+
    facet_grid(Year ~ type)+
    theme(legend.title = element_blank())+
    #scale_fill_viridis(discrete=TRUE,option = "")+
    #scale_x_continuous(breaks=seq(0,82,1))+
    ggtitle("TN sources")+
    ylab(bquote('TN x '* ~ 10^5*' (kg/yr)'))+xlab(paste(xlab))+
    theme(legend.title = element_blank(),axis.text.x = element_text(size=10),) 
  
}

graph1994(TN.1995,"1995")
graph3years(TN.yr,"")

#####EXport the 3 in one graph
t=Sys.Date()
png(filename=paste("Q:/Shared drives/Jordan Lake Modeling/WatershedModeling/Manuscript/Figures/Bar charts/TN 3basins_compared",t,".png",sep=""), 
    units="cm",width=16.5,height=14,res=300)

graph3years(TN.yr,"")
#graph3years_retvsnone(TN.total,"")

dev.off()

######ADD another 3 panels with 



total2=total %>%
  group_by(basin, Year) %>%
  summarise_each(funs(sum))
total2=as.data.frame(total2)
